Time-dependent bond-current functional theory for lattice Hamiltonians: 
fundamental theorem and application to electron transport 
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The cornerstone of time-dependent (TD) density functional theory (DFT), the Runge-Gross the- 
orem, proves a one-to-one correspondence between TD potentials and TD densities of continuum 
Hamiltonians. In all practical implementations, however, the basis set is discrete and the system 
is effectively described by a lattice Hamiltonian. We point out the difficulties of generalizing the 
Runge-Groos proof to the discrete case and thereby endorse the recently proposed TD bond-current 
functional theory (BCFT) as a viable alternative. TDBCFT is based on a one-to-one correspon- 
dence between TD Peierl's phases and TD bond-currents of lattice systems. We apply the TDBCFT 
formalism to electronic transport through a simple interacting device weakly coupled to two biased 
non-interacting leads. We employ Kohn-Sham Peierls phases which are discontinuous functions of 
the density, a crucial property to describe Goulomb blockade. As shown by explicit time propaga- 
tions, the discontinuity may prevent the biased system from ever reaching a steady state. 

PACS numbers: 31.15.ec, 73.23.Hk, 05.60.Gg 
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I. INTRODUCTION 

The central idea of timc-dcpcndent (TD) density func- 
tional theory (DFT) is to map a time-dependent and 
interacting many-particle system onto a time-dependent 
system of non-interacting particles moving in an effective 
Kohn-Sham (KS) potential, t;Ks(r,t), chosen such that 
the TD densities of the interacting and non-interacting 
systems are equal. From the knowledge of WKs(r,t) it 
is then possible to compute the TD density n(r, t), and 
hence the TD longitudinal current, in a one-particle man- 
ner. With this results as a starting point, an exact 
TDDFT formulation of quantum transport which can 
cope with both transient and steady-state regimes has 
been proposed [1, 2]. 

The theoretical foundation for the aforementioned 
mapping is laid down by the celebrated Runge-Gross 
theorem [3] whose essential statement is that the time 
evolution of two systems evolving from the same initial 
state 1 5*0) under the influence of two different TD poten- 
tials w(r,i) and v'{v,t), which are analytic in time and 
differ by more than a purely time-dependent, position- 
independent function C{t), leads to different TD densities 
ri(r,t) and n'{v,t). This theorem was developed further 
by van Leeuwen [4] who showed that for a many-body sys- 
tem with a given particle- particle interaction it;(|r — r'|) 
and moving in some TD potential w(r, t), the TD density 
can be reproduced in another many-body system with a 
different interaction w'dr — r'|) and moving in a TD po- 
tential u'(r,i) provided that the initial states |4'(0)) and 
|4''(0)) of the two systems yield the same density and lon- 
gitudinal current. Moreover, for a given initial state the 
potential w'(r, t) is unique up to a purely time-dependent 



function. Recently, Ruggenthaler and van Leeuwen ex- 
tended the validity of these results by relaxing the con- 
dition of Taylor-expandability in time [5] . 

The van Leeuwen theorem reduces to the Runge-Gross 
theorem for w = w' and |4'(0)) = |^''(0)), and estab- 
lishes the existence of a non-interacting KS system for 
w' = and |^''(0)) = |$ks) a Slater determinant. In 
this latter case the potential v'{r,t) reproducing a given 
density becomes the Kohn-Sham (KS) potential WKs(r, t) 
of TDDFT. In general, the KS potential is uniquely de- 
termined by the time-dependent density n{r,t), and the 
initial states |5'(0)) and |$ks)- If we further assume that 
the initial states |vE'(0)) and |$ks) arc ground states, then 
via the Hohenberg-Kohn [6] and Kohn-Sham [7] theorems 
of static density functional theory both initial states are 
uniquely determined by the initial density n(r, 0), and 
the KS potential ^;Ks(r,i) becomes a unique functional 
of the density alone. The TD density can then be calcu- 
lated from the TD KS equations [3] 
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where the sum in Eq. (2) is over the occupied orbitals in 
the KS Slater determinant. 

Generally, in practical implementations of the con- 
tinuum TD KS equations there are three independent 
sources of errors: (i) the use of only a finite number of 
basis functions, i.e., a de facto description of the problem 
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by a discrete or lattice Hamiltonian; (ii) the approxima- 
tion employed for the exchange-correlation (XC) poten- 
tial and (iii) inherent numerical inaccuracies. Even if we 
could use the exact KS potential (designed for contin- 
uum systems) and run the program on a machine with 
infinite accuracy, when using a finite basis set the den- 
sities of the KS system and of the original interacting 
system will, in general, differ. In fact, it would be useful 
to have a discrete version of TDDFT for ruling out the 
errors due to the incompleteness of the basis set. An even 
stronger motivation for constructing a discrete TDDFT 
is to have an alternative tool for the study of model sys- 
tems like, e.g., the Anderson model, the Hubbard model, 
etc. Recent works have shown that not all densities can 
be represented by lattice Hamiltonians [8]. The break- 
down of this u-reprcscntability is due to the fact that 
the rate of density change which can be supported on 
a given lattice is limited by the inverse of the hopping 
matrix elements [9, 10]. In the next section we point 
out an even more severe problem related to the difhcul- 
ties of generalizing the Runge-Gross proof to discrete or 
lattice Hamiltonians and re-understand why the contin- 
uum formulation works from the "discrete perspective" . 
To overcome these difficulties a TD bond-current func- 
tional theory (BCFT) has recently been proposed [13]; 
here the underlying idea is to use the so called Peierls 
phases (discrete version of the vector potential) instead 
of the scalar potential as the basic KS field. We will 
present the TDBCFT and demonstrate the analogous of 
the van Leeuwen theorem in discrete systems. In Sec- 
tion III we take advantage of the TDBCFT formulation 
to construct a TD framework for quantum transport in 
model Hamiltonians. As an application in Section IV 
we will revisit the Coulomb Blockade (CB) phenomenon 
in the light of a recent finding on the relevance of the 
derivative discontinuity of static DFT functionals [11] in 
the CB regime [12]. Conclusions and future perspectives 
are drawn in Section V. 



II. TIME-DEPENDENT DENSITY 
FUNCTIONAL THEORY ON A LATTICE 

A. Problems in generalizing Runge-Gross proof 

We write the Hamiltonian of interacting electrons in 
some orthonormal basis (pm{'<^) of orbitals localized at 
site m by expanding the field operators as tpa-(j) = 
J^m'^mafmiY), whcrc thc crcatiou (annihilation) oper- 
ators cl^^ {cm<y) satisfy the anticommutation relations 
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For electrons exposed to a TD on-site potential Vm{t), 
the Hamiltonian then reads 



where the non-interacting part consists of a static part 



and a time-dependent external potential 
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where we assume that v„i{t) = for all times t < 0. 
Hint denotes the operator describing the electron-electron 
interaction. Note that Vm{t) is coupled to thc density 
operator hm = So- ^lia^na at site to. 

In order to have a theoretical foundation for TDDFT 
on a lattice, one first has to prove the existence of a 
unique mapping between the expectation value Umit) of 
the density operator at site to and the on-site potential 
Vm{t)- The obvious idea is to adapt the original Runge- 
Gross proof [3] for continuum Hamiltonians to the case 
of lattice Hamiltonians. However, this strategy runs into 
problems as we will now demonstrate. 

Let l^'(t)) be the many-body state that solves thc TD 
Schrodinger equation 



^-\^it))^Hit)\^{t)) 
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with the initial condition |4'(t = 0)) = |^o)- Then the ex- 
pectation value 0{t) = {d{t)) = (*(t)|0(t)|*(t)) of any 
quantum mechanical operator 0{t) satisfies thc equation 
of motion 
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For a straightforward generalization of the Runge-Gross 
proof to lattice systems one needs to show that two wave- 
functions l^'(t)) and |^''(t)) which evolve under the influ- 
ence of two different potentials Vm (t) ^ (t) from the 
same initial state |\E'(0)) = |vl/'(0)) = |*o) always lead to 
different TD densities nm{t) ^ n'^it). Let us, for simplic- 
ity, consider the simplest case of two time-independent 
potentials Vm and v'^ that differ by more than an addi- 
tive constant. To show that they generate two different 
TD densities we calculate thc succesive time-derivatives 
in t = of rimit) and n'^(t) until the order (if any) at 
which they differ. The zero-th and first derivatives in 
t = are clearly the same. From Eq. (8) the second 
derivative of n,„ (t) and (t) in t = reads 



"-m(O) 



k 



Pkm H~ Pmk Tkm){vk (9) 



^'mW = ^{TmkPkni + P7nkTkm){v'k ^ "m) + A^, (10) 

k 

with p„ik = I]o-(*ol4cr Cmtrl^'o) th6 onc-particlc density 
matrix at t = and 
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cp(2)=0 



(p(3)=-l 

FIG. 1: The 4-site ring described in the main text with a 
single-particle eigenstate of energy zero. Any external po- 
tential Vm,{t) with vi{t) = V3{t) — leaves the density un- 
changed. 



focus on one-dimensional systems. Let m = 0, ±1, ±2, . . . 
be the label of a grid point Xm = rnAx of a one- 
dimensional grid with lattice spacing A^. Choosing the 
Tmn such that it gives the three-point discretization of 
the (continuum) kinetic energy 
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it is easy to show that the non-interacting Hamiltonian 
Ho{t) = Kq + V{t) (sec Eqs. (5) and (6)) approaches the 
continuum Hamiltonian 
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where we defined the field operators 
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Introducing the short-hand notation Kmk ~ T„ikPkm + 
PmkTkm for the kinetic energy density of the bond k — m, 
the difference Sri = n(0) — n'(0) is simply given by the 
formula below 
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where dv = v — v'. Like in the original Runge-Gross proof 
we now proceed by reductio ad ahsurdum. Assume that 
5n vanishes for a non-vanishing (and non-constant) 5v. 
Then, multiplying both sides of Eq. (12) by 5vm and 
summing over m we find 
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where we made use of the symmetry Kmk = Kkm ■ Due 
to the fact that the K^k do not have a definite sign Eq. 
(13) is not an absurdum. In fact, in a lattice theory we 
could have two different potentials that yield the same 
second derivative of the density in i = 0. We should then 
proceed further and try to prove that the third deriva- 
tives are different. One immediately realizes, however, 
that this way does not lead anywhere. There exist several 
counterexamples to the existence of a one-to-one corre- 
spondence between densities and potentials in a lattice 
Hamiltonian. Consider, for example, a one-dimensional 
ring with four sites, zero on-site energies Tmm = and 
nearest-neighbor hopping T12 = T23 + T34 = T41 = T, 
see Fig. 1. The single-particle state with amplitudes 
(p{l) = 1, (f{3) = — 1 and (p{2) = (p{4) = is an eigen- 
state with energy zero. Any external potential Vm (t) with 
vi{t) = V3{t) = leaves the density unchanged. 

In the next Section we show how to overcome these dif- 
ficulties by considering the Peierls phases instead of the 
on-site potentials as the basic KS fields. Before, how- 
ever, we find particularly instructive to re-understand 
why the continuum formulation works starting from the 
"discrete" perspective. To limit the complications we will 



and the scalar potential v(xm,t) = Vm{t). Inserting the 
T,nn of Eq. (14) into Eq. (13) and taking the limit -> 
we find 
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with n{x) = limA^-j-o R-e 
sity. For all situations of physical interest the integrand 
may vanish at most in a set of zero measure and is oth- 
erwise larger than zero. Thus, in the continuum case the 
above equation is an absurdum, in agreement with the 
Rungc-Gross theorem. 



B. TDCDFT on a lattice 

The problems in extending the Runge-Gross theorem 
to lattice systems sketched in the previous Section can be 
circumvented by using TD Peierls phases instead of TD 
on-site potentials [13]. For arbitrary TD electromagnetic 
fields described by a scalar and vector potential v{r,t) 
and A(r, t) one can always perform a gauge transforma- 
tion such that v{r,t) — >■ v{r) and A{r,t) — ^ A(r,t), so 
that the time-dependence resides only in the vector po- 
tential. In particular, for those situations when there 
is only a TD potential w(r,<), we could always gauge v 
away and work with a TD vector potential A(r, t). In the 
localized basis we choose, therefore, the time-dependent 
Hamiltonian to have the form 



where 



H{t) = K{t) + H,,,t, (18) 
^(*) = E E Tmne'-'-^'^cl^Cr,,, (19) 
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and the time-dependence only resides in the Peierls 
phases 7mn(0 = —lnm{t)- In a grid basis with grid 
points r„i, the on-site matrix elements Tjnm ~ v{r„i) 
contain information on the on-site electrostatic potential 
v{r) while the phases jmn{t) = \ J^" dl ■ A(r, t) describe 
the effects of an external TD vector potential A(r,t). 

From Eq. (8) the density at site m satisfies the equa- 
tion of motion 

-^n„i(t) = ^ Jrnn{t) + i {| [-f^int ) (20) 
n 

where we have defined the bond-current operator Jm„ (t) 
according to 

J™«(i)-|^(T„„e'T-(*)cLc„.-H.c.) . (21) 

(T 

In turn, the equation of motion for Jmn{t) reads 

^Jmnit) = i^™„(i)^7m„(i) + Fm„ (0 (22) 

with the "kinetic energy" density operator 

Kmnit) = J2 {T-mne'-'-^'hl^Cr,^ + H.C.) (23) 

and 

{t)^t[H{t),J„uM- (24) 

We now ask the question whether the bond-currents 
Jmn{t) along those bonds with T^n can be repro- 
duced by another Hamiltonian H'{t) with a different in- 
teraction and a different TD electromagnetic field, i.e., 

H'{t) = k'{t)+HU (25) 

with 

k\t) = Tmne'^'—^'^l^Cn. • (26) 

Note that the matrix elements Tmn in Eq- (26) are the 
same as in Eq. (19). Let us define the initial configuration 
of the primed system as the couple {|*I'o), 7'(0)}, that 
consists of the initial state and the initial value of the 
Peierls phases. For the current J^„{t) of the primed 
system to be the same as the current Jm„ (t) at the initial 
time t = 0, the initial configuration must be compatible^ 
i.e., it must fulfill 

(*ol| E (7^™"e*^-WcL5«. - H.C.) 1*^,) = J™„(0). 

(27) 

Assuming the existence of at least one compatible config- 
uration we demonstrate below that there exist TD Peierls 
phases {Y{t)} for which Jmn{t) = J'mni^) for bonds 
with T'mn(O) 7^ 0, and that these phases are unique. 



The bond current densities Jmn{t) and J^„j(i) of the 
two systems are the same provided that 

Am„(0^7™„(0 =^™4i)^7m«(0+^™«(0-Km(0 • 

(28) 

To proceed further we adopt the same strategy of Vignale 
[14] in the context of TD current DFT, and expand all 
quantities in Eq. (28) in a Taylor series around t = 0, 
e.g., 

OO 

K,„{t)=Y,K'^lt' . (29) 

1=0 

Inserting the Taylor expansions into Eq. (28) and equat- 
ing the coefficients with power / then gives 

K'^'r! {I + ih'^:'' = - E(fc + m'^u'hii^'' 

fe=0 

+E(fc + i)^il;SL+^^ + ^^i^I-^^^2- (30) 

fe=0 

The quantities K'^^{t) and F^^{t) depend explicitly on 
the phases {j'{t)} and therefore their Z-th derivative de- 
pend on all the fc-th derivatives {^'C^)} with k < I. Thus, 
under the mild condition that the initial compatible con- 
figuration is chosen in such a way that 

if:„„(0) = K'^°2 + , (31) 

Eq. (30) constitute a system of recursive relations to 
determine all Taylor coefficients of 7m„(i) for all bonds 
with Tmn 7^ 0. 

A direct corollary of this result is that if the inter- 
action Hamiltonian = Hmt and the initial config- 
uration {|*o)'7'(0)} = {l*o),7(0)}, then the Peierls 
phases 7'(t) = 7(i), i.e., for any fixed initial configu- 
ration there is a one-to-one correspondence between the 
bond-currents and the Peierls phases. Consequently, we 
can think of the TD many-body state |^'(i))j and hence 
of all TD expectation values, as a functional of the bond- 
currents: l^'(t)) = |^'[J](t)). Another important corol- 
lary is that the bond-currents of a system with interac- 
tion iJint can be reproduced in a system of noninteracting 
electrons, Hl^^ = 0, which we call the KS system. These 
two corollaries lay down the foundations of a TD bond- 
current functional theory (BCFT) for discrete (or lattice) 
Hamiltonians, where the basic variable is the bond cur- 
rents J{t) and the basic KS field is the Peierls phases 
Jit). 

Once the equality of the bond currents is established, 
it follows from Eq. (20) that if 

[Hint.flrn] = [H-^t^fljn] = (32) 

the time evolution of the densities nm{t) and n'^{t) is 
also the same, provided that they are the same at the 
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initial time, nm(0) = n'^ (O). It is worth noting that for 
Hubbard-hke interactions, Eq. (32) is satisfied. 

Finally we would like to emphasize that in the above 
proof the possibility of reproducing the bond-currents in 
a system with a different interaction heavily relies on the 
convergence of the Taylor series, which has here simply 
been assumed [15]. An alternative proof based on the 
Picard-Lindelof theorem for non-linear differential equa- 
tions has been recently proposed by Tokatly [16]; this 
proof has the advantage of avoiding the Taylor expan- 
sions but it remains unpredictive on the maximum time 
until which a solution exist. Another merit of Tokatly's 
stategy is that TDBCFT can be generalized to primed 
systems with different hopping integrals ^ Tmn- 

III. MODEL HAMILTONIAN FOR 
TIME-DEPENDENT TRANSPORT 

Tight-binding based model systems of interacting elec- 
trons are often used to study transport through nanos- 
tructures such as quantum dots. Although most of these 
studies are restricted to the steady-state regime, more 
recently there has been increasing activity to describe 
the time evolution towards the steady state as the sys- 
tem is driven out of equilibrium by applying a bias in 
the leads. These studies use a range of methods such 
as, e.g., TDDFT [2, 17-20], generalized master equations 
[21], many-body perturbation theory [22-24], the time- 
dependent density-matrix renormalization group [25-27], 
a quantum trajectory approach [28], or real-time path in- 
tegral [29] and Monte Carlo approaches [30]. 

Here we will apply the TDBCFT formalism developed 
in the previous Section to describe TD quantum trans- 
port through a one-dimensional Anderson-like model sys- 
tem [31]. We will consider a two-terminal setup where 
two semi-infinite non-interacting leads are coupled to an 
interacting single-level impurity. The total Hamiltonian 
of the system reads 

Hit)^Hc{t) + Hi^it) + Hn{t) + HT. (33) 
We model the Hamiltonian for the impurity as 

Hcit) = ^cit)dld, + Hint, (34) 

with a Hubbard-like electron-electron interaction 

Hint = Ud\d^dldi. (35) 

The Hamiltonian for the semi-infinite lead a = L, R is 
taken of the form 

oo 

Z— 1 (T 

OO 

-EE (^"C,Lc.+i.o + H.c.) , (36) 

i—1 (7 



and the tunneling Hamiltonian 

Ht = -Y. E (l^ink^Cl.a + h.c.) (37) 
(T a=L,R 

connects the impurity to both left and right leads. 

We note that the Hubbard interaction (35) satisfies 
Eq. (32) since it commutes with the on-site density fi^ = 
dl-da in the impurity as well as with the on-site density 

'^maa — ^rnaa^^^^cv lead Oi, i.C, 

[Hint,n„] = [Hint^nmaa] = . (38) 

In our model the time-dependence appears exclusively 
through a TD on-site potential, that is the bias Wait) 
in lead a and the on-site energy ec{t) at the impurity. 
Therefore, the gauge transformation Cmaa — ^ Cmo-aC*^"'*^ 
with Pa{t) = J^dt'Wa{t') and d^ d^e*'^c(t) ^ith 

l^cit) = /o di'ec(*') leads to Peierls phases of the par- 
ticular form ^rnanait) = for any pair of sites m and n 
in lead a and to 

7L(i)= f dt' (WUf) - ec{t')) (39) 
Jo 

and 

7R(t)= f dt'{ec{t')-Wi^{t')) (40) 
Jo 

for the Peierls phases of the bonds connecting the left 
and right leads to the impurity. Under this gauge trans- 
formation the TD on-site potential is gauged away and 
the Hamiltonian (33) is transformed into an Hamiltonian 
of the form in Eqs. (18-19) 

H{t)^Hc + Hi^ + H^ + HT{t). (41) 

The transformed impurity Hamiltonian Hq = Hi^t and 
lead Hamiltonians 

_^ oo oo 

-f^a = E E E (^«ctaCi+laa + H.C.) 

(42) 

become independent of time whereas the transformed 
tunneling Hamiltonian 

Mt) = -Yl E (^inke^^°^*^4ciaa+H.c.) (43) 
carries all the time dependence. 

A. Local Peierls phase approximation and the 
ABALDA functional 

The TDBCFT presented in Section II B guarantees 
that the bond-currents Jmn of the interacting system can 
be reproduced in a non-interacting KS system with TD 
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FIG. 2: Schematic representation of the model system. The 
biased system (top) with Hubbard-type interaction on the 
dot is gauge-equivalent to the system (middle) with time- 
dependent Peierls phases 7i(t) and 7ii(t) (see Eqs. (39) 
and (40)). In the local Peierls phase approximation, the 
Kohn-Sham system (bottom) is assumed to have nonvanish- 
ing Peierls phases (Eqs. (44) and (45)) only at the same links 
as the interaction system (middle). 



Peierls phases ^^ni^)- In the following we make a lo- 
cal approximation for the Jmni^)- set them to zero 
except for those bonds where the external Peierls phases 
lmn{t) are non-zero, i.e., the bonds connecting the impu- 
rity to the leads, see Fig. 2. Within this approximation, 
we write the only two KS Peierls phases similarly to Eqs. 
(39-40) 



di'(M^L(i')-«Ks(t')): 



and 



7rw= / di'(^,Ks(t')-vi^L(0): 



(44) 



(45) 



where ec (t) is replaced by the on-site KS potential f ks (t) 
that we write as the sum of the external on-site poten- 
tial ec{t), the Harteee potential ^U'n{t) (with n(<) the 
TD density at the impurity) and the exchange-correlation 
(XC) potential Wxc(i) 



VKs{t) = Eo (0 + 7;Uno{t) + Uxc(i) ■ 



(46) 



Note that by a gauge transformation we could go back to 
the original gauge in which there is an on-site potential 
instead of the Peierls phases. In the original gauge the 
KS system is described by the Hamiltonian 

H{t) = ff^S(t) + ^L(t) + Hnit) + Ht, (47) 

where the lead and the tunneling Hamiltonians, Ha{t) 
and Ht, are unchanged compared to the interacting 
Hamiltonian (i.e., they are given by Eqs. (36) and (37), 
respectively) and the KS Hamiltonian for the impurity is 



KS 



(48) 



We still need to specify the approximation used for 
the XC potential i'xc(i)- In this work we use an adia- 
batic version of the local density approximation (LDA) 
for the static, non-uniform Hubbard model. The con- 
struction of this functional [32] follows a similar strategy 
as the one used to construct the usual LDA based on 
the model of the uniform electron gas. The crucial dif- 
ference is that the underlying model is the uniform Hub- 
bard model whose exact solution can be constructed via 
the Bethe ansatz. Just as the XC energy of the uniform 
electron gas serves as input for the usual LDA, the exact 
XC energy per site of the uniform Hubbard model then 
serves as input for the Bethe-ansatz LDA (BALDA) used 
for non- uniform Hubbard models [32, 33]. In the context 
of TDDFT, an adiabatic version of this functional (adia- 
batic BALDA, ABALDA) local in both space and time, 
e.g., Vxc["-](*,i) ~ Wxc("-i(i)), has been suggested by Ver- 
dozzi [10]. The form of the BALDA XC potential derives 
from a parametrization of the XC energy per particle of 
the uniform Hubbard model. This parametrization de- 
viates somewhat from the exact, numerical XC energy 
of the Hubbard model [34, 35], especially for weak in- 
teractions, but here we are not concerned about these 
differences. The crucial property for our purposes is the 
existence of a derivative discontinuity at half filling [36] 
(see discussion below). 

The original BALDA has been proposed for a system 
which not only has the same on-site energy but also the 
same interaction for all sites and also the hopping con- 
necting neighbouring sites is the same throughout the 
chain. In Ref. 33 the BALDA functional has already 
been used successfully for site- dependent interactions. It 
should be noted that this is a deviation from the usual 
LDA philosophy where the XC energy of the uniform gas 
is used for a non-uniform system but the interaction is 
unchanged. Here we also take the same approach and use 
the BALDA for a system with interactions only at one 
site. In addition, we are interested in situations when 
the interacting region is weakly connected to two non- 
interacting leads, i.e., the hopping V within the leads 
may be different from the coupling Viink of the leads to the 
interacting region. For this situation a generalization of 
the original BALDA functional has been suggested [12]. 
The explicit form of this modified BALDA XC potential 
reads 

v^c[n] = 9(1 - n)vi<^[n] - e{n - l)ui<)[2 - n], (49) 
where 



4c^M = -^C^«-2yn„k cos(^ 



\ / nn 

— cos 

2 J U 



The parameter ^ is determined by the equation 

Jq{x)Ji{x) 



2Q f-oa 

— sin(7r/C) = 4 dx — 
TT Jo a;[l+exp(t/x/(2Mi„k))] 



(50) 
(51) 



where Ji=o.i{x) are Bessel functions. The BALDA XC 
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FIG. 3: Smoothened BALDA exchange-correlation potential 
(Eq. (53)) as function of the density for various values of 
U /Viiny. with smoothing parameter a — 10^'*. 



potential has a discontinuous jump at half-filling [36]: 

?;xc(n = l+)-Uxc(ri = 1") = ?/-4ViinkCos . (52) 

Here it is interesting to note that in the limit of very 
weak coupling (Vlink 0) the discontinuity just reduces 
to U which is the charging energy required to put a sec- 
ond electron on the interacting impurity if it is already 
occupied by one electron, i.e., in this limit the functional 
becomes exact. Although the exact XC potential for the 
Hubbard model is certainly discontinuous, in our situ- 
ation where a single interacting impurity is coupled to 
two non-interacting leads, the coupling to the leads in- 
troduces some broadening in the levels of the isolated 
impurity and it is reasonable to expect that this broad- 
ening leads to a smoothening of the discontinuity. We 
therefore introduce a smoothening in our model and use 
the following XC potential 

«xcN = /(«)4c'N - (1 - /(«))4c'[2 - n\ 



(53) 



where 



l)la) 



1 



(54) 



exp((n 

and a is a smoothening parameter. In Fig. 3 we show 
the BALDA XC potential for different values of [//Mink 
using the value a = 10~* for the smoothening parameter. 

It has recently been shown [12] that this discontinuity 
is crucial to describe Coulomb blockade. Moreover, in a 
time-dependent picture and in the parameter regime of 
Coulomb blockade it also prevents the biased system to 
reach a steady state. More details on these findings will 
be discussed below. 



B. Time propagation with embedding 

By mapping the interacting problem onto a non- 
interacting one we have already achieved a considerable 



simplification. However, we are still dealing with an in- 
finitely extended system which has to be treated numer- 
ically. This can be achieved by an embedding technique 
which maps the infinite system exactly onto a tractable, 
finite problem. 

In the localized site basis the KS Hamiltonian (48) has 
the matrix structure 



(55) 



It is possible to derive the equation of motion for the k- 
th KS single-particle orbital projected onto the central 
region, 'ipk,c{t)^ which reads 
















H = 












V 







V'fc,cW= / dt'Se,n(t,f')'/'fe.c(i') 

Jo 

+ ^Hc„gac.(i,0)i^fc,„(0),(56) 



where 



^cni{t,t') — ^ llcaSaa{t,t')'H.aC 



(57) 



a=L,R 



is the retarded embedding self energy. Here, the retarded 
Green's function gaa{t,t') of the isolated lead a satisfies 
the equation of motion 



^^-H„„(t) 



(58) 



with boundary conditions Saait,t^) ~ and 
gfjQ,(t,t~) = —i. The TD density at the impurity can 
be calculated from the KS orbitals as in Eq. (2), i.e.. 



(59) 



where the sum runs over the occupied KS states. In the 
present work we use the algorithm described in Ref. 2 
to propagate Eq. (56) starting from the self-consistent 
BALDA ground state of the coupled lead-dot-lead sys- 
tem. 



IV. THE DERIVATIVE DISCONTINUITY AND 
ITS CONNECTION TO COULOMB BLOCKADE 

A. Ground state and non-equilibrium steady state 

We will use the KS Hamiltonian in the BALDA ap- 
proximation to describe transport for our model of a sin- 
gle interacting impurity connected to two non-interacting 
tight-binding leads. This model has been studied exten- 
sively in the literature, and the results from other meth- 
ods can be used for a validation of our KS treatment. 
In the present Section we first assess the quality of our 



approximate BALDA functional in the ground state by 
comparing to recent Quantum Monte Carlo (QMC) re- 
sults for the same model system. Then we look into the 
nonequilibrium steady state for the system exposed to a 
DC bias in the leads, assuming that such a steady state 
exists. 

Using the techniques of non-equilibrium Green's func- 
tions one can derive a self-consistency condition for the 
ground or steady-state electron density n°° at the impu- 
rity which reads 



= 2 E 



a=L,R 



^r{u-W^)\G{u)\' . (60) 



Here, Wa {a = L, R) is the constant bias applied in lead 
a (of course, the ground state density is obtained for 
VFl = Wr = 0) and 



G(w) 



uj - VKS [n 



(61) 



=L3 



is the retarded Green's function in frequency space at the 
impurity. T,a{i-u) is the Fourier transform of the a-lead 
contribution to the embedding self energy (57) which, 
for the present case of non-interacting tight-binding leads 
with vanishing onsite energies, = 0, is given explicitly 
by 




for bJa > 2Vq 

for uja < ~2Va (62) 

for \uJa\ < 2Va 



with real and imaginary parts Aa{uj) and Ta{uj), re- 
spectively. Finally, the total width function is r(a;) ~ 
J2a^a{i^) and is the Fermi energy of the contacted 
system in the ground state. 

We have first solved Eq. (60) for the ground state den- 
sity no := n°° for Wl = = and compared to 
recent QMC calculations presented by Wang et al. in 
Rcf. 37. Here we focus on the dependence of uq on 
the on-site energy ec at the impurity. We set e_F = 
(half filling of the leads) and for the Hubbard interac- 
tion we use U = 0.105,0.21,0.42,0.84 (ah energies are 
given in units of the equal hopping in left and right leads 
Vl — Vr = V). For the hopping between leads and im- 
purity we take a weak Mink = 0.1803. In order to achieve 
a meaningful comparison to the QMC data, the value of 
our Viink is a factor of l/\/2 smaller than the one used 
in Ref. 37 since Wang et al. considered an impurity cou- 
pled to a single lead only [38]. For the smoothening pa- 
rameter we choose the value a = 10^^. Fig. 4 shows the 
ground state density as function of on-site energy ec- We 
see that the QMC and BALDA results agree surprisingly 
well, especially for weak interactions. For strong inter- 
actions the BALDA develops a clear Coulomb blockade 




FIG. 4: Comparison of BALDA and Quantum Monte Carlo 
(QMC) ground state densities at the impurity as function of 
onsite energy Eq for different values of the interaction. The 
coupling between leads and impurity is Viink = 0.1803. All 
energies are given in units of the hopping V in the leads. The 
QMC results were extracted from Ref. 37. 



step, i.e., the density hardly changes over a significant 
range of onsite energies. Although the step in BALDA 
extends over a smaller range of onsite energies than in 
QMC, the agreement is still quite reasonable. Already at 
this stage it is clear that the feature which gives rise to 
the Coulomb blockade step is the derivative discontinuity 
built into the BALDA functional. 

Now we turn our attention to the biased, non- 
equilibrium situation. If we assume that the system 
evolves toward a steady state in the long time limit then 
the steady-state value of the density at the impurity is 
given by the self-consistent solution of Eq. (60). We solve 
the latter for different values of the bias VFl applied in 
the left lead. The system parameters are ec = 2, U = 2, 
ep = 1-5 and a = 10~^. 

In the left panel of Fig. 5 we show the steady state 
density as function of the applied bias for different val- 
ues of the hopping Viink between leads and impurity. We 
clearly see a plateau in the steady-state density at unity, 
i.e., when the impurity is occupied by exactly one elec- 
tron. This plateau becomes wider and more step-like as 
Mink decreases which corresponds to the usual picture of 
Coulomb blockade: the dot can only be occupied by zero, 
one, or two electrons, and the energy cost for double oc- 
cupancy is given by the Hubbard interaction U . It is im- 
portant to emphasize here that in the plateau region it is 
crucial that we used an XC potential with a smoothened 
discontinuity. If one uses the truly discontinuous XC po- 
tential, the steady-state condition actually does not have 
any solution at all in this parameter region, indicating 
that the initial steady-state assumption was not justi- 
fied. This point will be further investigated in the next 
Section. 

In the right panel of Fig. 5 we show the steady-state 
densities for the same parameter range if one uses the 
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FIG. 5: Steady state density as function of the applied bias 
Wl for different values of Vunk- Left panel: BALDA steady 
state densities, right panel: Hartree steady state densities. 

Hartrcc approximation, i.e., when Wxc in Eq. (46) is set 
to zero. One clearly sees that in this case the step in 
the steady-state density is completely absent. Instead, 
for small Vunk the steady-state density rises almost lin- 
early within a certain bias range. For larger values of 
Mink, the Hartree and BALDA steady-state densities are 
qualitatively quite similar. 

B. Time-dependent transport 

In the previous Section we assumed that under appli- 
cation of a DC bias the system evolves towards a steady 
state for which we then solved a self-consistency condi- 
tion to obtain the steady-state density. In the present 
Section we will show that for certain parameters the 
steady-state assumption is actually not always justified: 
when perturbing the system, which is initially in its 
ground state, by switching on a DC bias the time evolu- 
tion does not necessarily lead to a steady state but rather 
to an oscillatory state of dynamical density oscillations. 

In Fig. 6 we show the time evolution of the quantum- 
dot density from its intial, ground-state value when, at 
t = 0, a. DC bias is suddenly switched on in the left 
lead, i.e., W^it) = 6i(i)WL with Wl = 0.5,1.3,1.6,1.9. 
The parameters are the same as those used for Fig. 5 
with Viink = 0.3. According to the steady-state picture 
of Fig. 5, the first value Wl = 0.5 corresponds to the 
rising flank of the density as function of bias while the 
other biases correspond to the plateau. The TD density 
corresponding to bias = 0.5 can be seen to evolve 
smoothly towards its steady state value. However, look- 
ing closely at the evolution of the density for the other 
bias values (inset of Fig. 6) one can see that for these 
cases the system does not evolve towards a steady state 
but rather approaches a dynamical state of periodically 
oscillating density. Furthermore, the density oscillates 
around unity, i.e., single occupancy of the impurity. The 
amplitude of these oscillations is rather small, of the or- 



— Wl = 0.5 

— 1.3 - 
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FIG. 6: TD density in BALDA for different values of the bias 
Wl for Viink ~ 0.3. Inset: Magnification of the TD density at 
large times. The straight line at unity is a guide for the eye. 

der of 5 X 10^'^. We also note that the larger the value 
of the bias (in the range of the plateau of Fig. 5), the 
larger the fraction of time for which the density exceeds 
the critical value of unity. 

The qualitative difference of the response of the sys- 
tem for biases within and outside the step region of 
Fig. 5 becomes quite apparent when looking at the time- 
dependent KS potentials (left panels of Fig. 7) and cur- 
rents through the impurity (right panels of the same fig- 
ure). While for the bias outside the step region (Wl = 
0.5) both quantities behave smoothly, for the other bi- 
ases the KS potentials show rapid variations which are 
due to the discontinuity at n = 1: as long as the time- 
dependent density stays below or above unity, the KS 
potential changes smoothly in time. However, when the 
density crosses unity, the KS potential jumps discontinu- 
ously (or very rapidly for the smoothened KS potentials 
used here) . The rapid variation of the KS potential al- 
lows us to understand why for these cases a steady state 
is not achieved: consider the situation when the time- 
dependent density crosses the critical value of unity from 
below. At the instant td when the jump in the KS poten- 
tial occurs, the rate of change of the density is positive, 
n{tci) > 0, and immediately after td the density contin- 
ues to grow because the maximum rate of change of the 
density is limited by the inertia of the electrons. How- 
ever, for times right after td the significant increase of 
the KS potential tends to push down the density, i.e., 
n{t) < 0, and therefore n{t) will cross the critical value 
of unity from above at some time tc2 > td ■ At that time 
n{tc2) < 0, the KS potential will suddenly be lowered 
and electrons will be attracted by the quantum dot, i.e., 
n{t) > 0. Therefore, the density will eventually crosses 
unity at some time tc3 > ic2 from below and the cycle 
will start over. 

The above results have been obtained with an XC po- 
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FIG. 7: Time-dependent KS potentials and currents at the 
impurity for different values of the bias Wl for Viink = 0.3. 
Left panels: BALDA KS potentials, right panels: TD currents 
at the impurity site. 

tential with a smootheiicd discontinuity for which, in 
principle, the condition for the steady-state density dis- 
cussed in the previous Section does have a solution. This 
solution is not reached by the time propagation because 
the rate of change of the density is limited by T^ink [10]. 
The rate of change of the KS potential when the density 
crosses the (smoothened) discontinuity is instead limited 
by the smoothening parameter a~^. Therefore the steady 
state should not be attained whenever the rate of change 
of the KS potential is larger than the largest possible rate 
of change of the density. 

The time-dependent currents in the Coulomb blockade 
regime (lower three right panels of Fig. 7) show sawtooth- 
like oscillations at the impurity site. By the equation of 
motion for the currents, this is consistent with the fact 
that the KS potentials form a train of almost rectangular 
potential steps. In contrast, the current away from the 
Coulomb blockade regime (upper right panel of Fig. 7) 
evolves smoothly towards its steady state value. This 
is qualitatively similar to the behaviour which is found 
for time-dependent Hartree calculations (in any parame- 
ter regime). Finally, we would like to point out that the 
oscillations just described are entirely due to electron cor- 
relations and are therefore different from oscillations oc- 
curing due to the presence of single-particle bound states 
[39-41]. 

V. CONCLUSIONS 

In the present work we have described the difficulties in 
generalizing the fundamental Runge-Gross proof of time- 
dependent density functional theory to systems described 
by a lattice Hamiltonian. To overcome these difficulties 



we have employed the time-dependent bond current in- 
stead of the density as basic variable. We were then 
able to prove a one-to-one correspondence between time- 
dependent bond currents and Peierls phases describing 
the electromagnetic field on the lattice and therefore pro- 
posed a time-dependent bond current functional theory 
as the proper extension of TDDFT to lattice systems. 

In a second part we have described a time-dependent 
approach to electron transport in model systems using 
this TDBCFT. We proposed a local Peierls phase approx- 
imation assuming that the Peierls phases of the nonin- 
teracting KS system are nonvanishing only for the same 
links as for the interacting case. A simple model of an in- 
teracting impurity connected to two noninteracting leads 
has been studied employing a functional which exhibits 
an explicit derivative discontinuity leading to a time- 
dependent KS potential which jumps discontinuously as 
the particle number on the impurity crosses an integer 
[42]. It has been demonstrated that this discontinuity 
is (i) crucial to describe Coulomb blockade and (ii) pre- 
vents the system from reaching a steady state by time- 
evolution out of the initial ground state upon application 
of a bias in the leads. Instead we find that after some 
transient time the system reaches a dynamical state of 
undamped density oscillations. In this picture, Coulomb 
blockade manifests itself as a dynamic phenomenon of 
sequentially charging and discharging the impurity. 

Here we demonstrated the crucial role played by the 
discontinuity in the XC potential to properly describe 
Coulomb blockade, a physical phenomenon which is due 
to electronic correlations, for a simple model system. 
This model is not easily generalized to more complicated 
situations such as, e.g., more orbitals or degrees of free- 
dom per site or more complicated lattices. As for the lat- 
ter case, there have recently been activities to construct 
the XC energy per particle which also display a deriva- 
tive discontinuity for more complicated model systems 
such as a graphene-like hexagonal Hubbard lattice [43] 
or the Hubbard model in three dimensions [44]. How- 
ever, the challenge remains to construct practical and 
universal approximations with a derivative discontinuity 
not only for lattice systems but also for systems described 
by continuum Hamiltonians. 
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